rm(list=ls())
cat("\014") 

# install.packages(c('ordinal','doRNG','foreach','doMC'))

require(ordinal)
require(doRNG)
require(foreach)
require(doMC)

#################
### FUNCTIONS ###
#################

inv_logit=function(logit) {
  exp(logit)/(1+exp(logit))
}

runSIM=function(proposal,window,data,ord_beta,nonresponse_beta,N_responses,range=data.frame(cfscore=seq(-1.626,1.653,by=.01)),sim_no=1) {
  temp=data
  # gen population parameters
  b01=proposal+(window/2)
  b02=proposal-(window/2)
  b1=ord_beta
  # log-odds
  logodds1=b01+b1*temp$cfscore
  logodds2=b02+b1*temp$cfscore
  prob_2to3=inv_logit(logodds1)
  prob_3=inv_logit(logodds2)
  # remaining probs
  prob_1=1-prob_2to3
  prob_2=prob_2to3-prob_3
  # simulate DV
  for (i in 1:nrow(temp)) {temp$y[i]=sample(1:3,size = 1,prob=c(prob_1[i], prob_2[i], prob_3[i]))}
  temp$y=as.factor(temp$y)
  # proposal estimate with full response
  m=clm(y~cfscore,data=temp,family='probit')
  range$pr.prop=predict(m,range,type='prob')$fit[,'2']
  est_proposal=range$cfscore[range$pr.prop==max(range$pr.prop)]
  # proposal estimate with nonresponse
  obs.temp=temp[sample(nrow(temp),N_responses,prob=inv_logit(nonresponse_beta*temp$cfscore*(1-temp$firm))),]
  m_obs=clm(y~cfscore,data=obs.temp,family='probit')
  range$pr.prop.obs=predict(m_obs,range,type='prob')$fit[,'2']
  est_proposal_obs=range$cfscore[range$pr.prop.obs==max(range$pr.prop.obs)]
  # return simulation results
  return(data.frame(
    'Proposal'=proposal,
    'Window Size'=window,
    'Estimated Proposal'=est_proposal,
    'Estimated Proposal (Observed)'=est_proposal_obs,
    'Bias'=proposal-est_proposal_obs,
    'No. of Commenters'=N_responses,
    'Beta for CFscore'=ord_beta,
    'Beta for CFscore on Nonresponse'=nonresponse_beta,
    'Sim. Number'=sim_no))
}

#################
## SIMULATIONS ##
#################

# Iterate over Various Settings 
Ns=seq(100,500,by=100) # Number of Comments
Props=seq(-1,1,length.out=5) # Lib-Conservatism of Proposal
Windows=seq(0.5,1,length.out=4) # Size of true ''Just Right'' Response Window
Betas=seq(-0.5,-1,length.out=5) # Magnitude of Relationship between CFscore and Comment Content
Betas_NR=seq(0.25,1,length.out=4) # Magnitude of Relationship between CFscore and Non-Response

load('sim-pop.Rda')

sims=500
seed=1990

registerDoMC(cores=detectCores()-1)

settings=expand.grid(Ns,Props,Windows,Betas,Betas_NR,KEEP.OUT.ATTRS=F)
names(settings)=c('Ns','Props','Windows','Betas','Betas_NR')


# generate bias figures under various assumptions
t0=proc.time()
simulations=foreach(x=1:nrow(settings),.combine=rbind,.options.RNG=seed) %dorng% {
  results=data.frame()
  for (t in 1:sims) {
    result=runSIM(proposal=settings$Props[x],
                  window=settings$Windows[x],
                  data=pop,
                  ord_beta=settings$Betas[x],
                  nonresponse_beta=settings$Betas_NR[x],
                  N_responses=settings$Ns[x],
                  sim_no=t)
    results=rbind(results,result)
  }
  return(data.frame(
    'Proposal'=settings$Props[x],
    'Window'=settings$Windows[x],
    'Commenters'=settings$Ns[x],
    'Beta'=settings$Betas[x],
    'Beta_NR'=settings$Betas_NR[x],
    'Estimated_Proposal'=mean(results$Estimated.Proposal),
    'Estimated_Proposal_Observed'=mean(results$Estimated.Proposal..Observed.),
    'Bias'=mean(settings$Props[x]-results$Estimated.Proposal..Observed.),
    'SE_Bias'=sd(settings$Props[x]-results$Estimated.Proposal..Observed.)))
}
proc.time()-t0

# save(simulations,file='sim-results.Rda')